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ABSTRACT 

Context. Neutron stars with strong toroidal magnetic fields are often produced in nature. We show that isentropic neutron stars with 
purely toroidal magnetic fields are unstable against the interchange, Parker and/or Taylor instabilities irrespective of the toroidal 
magnetic field configurations. 

Aims. The aim of this paper is to clarify the stabilities of neutron stars with strong toroidal magnetic fields against non-axisymmetric 
perturbation. The motivation comes from the fact that super magnetized neutron stars of ~ lO'^G, magnetars, and magnetized proto- 
neutron stars bom after the magnetically-driven supemovae are likely to have such strong toroidal magnetic fields. 
Methods. Long-term, three-dimensional general relativistic magneto-hydrodynamic simulations are performed, preparing isentropic 
neutron stars with toroidal magnetic fields in equilibrium as initial conditions. To explore the effects of rotations on the stability, 
simulations are done for both non-rotating and rigidly rotating models. 

Results. We find the emergence of the Parker and/or Tayler instabilities in both the non-rotating and rotating models. For both non- 
rotating and rotating models, the Parker instability is the primary instability as predicted by the local linear perturbation analysis. 
The interchange instability also appears in the rotating models. It is found that rapid rotation is not enough to suppress the Parker 
instability, and this finding does not agree with the perturbation analysis. The reason for this is that rigidly and rapidly rotating stars 
are marginally stable, and hence, in the presence of stellar pulsations by which the rotational profile is deformed, unstable regions 
with negative gradient of angular momentum profile is developed. After the onset of the instabilities, a turbulence is excited. Contrary 
to the axisymmetric case, the magnetic fields never reach an equilibrium state after the development of the turbulence. 
Conclusions. Isentropic neutron stars with strong toroidal magnetic fields are likely to be always unstable against the Parker instability. 
A turbulence motion is induced and maintained for a long time. This conclusion is different from that in axisymmetric simulations 
and suggests that three-dimensional simulation is indispensable for exploring the formation of magnetars or prominence activities of 
magnetars such as giant flares. 
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1. Introduction are associated with proton cyclotron lines, the magnetic field 

l/^ . strength is estimated to give B ~ 10'^ G. 

' For about a dozen accreting X-ray pulsars in binary systems, 

O There are a lot of observational evidences that suggest the pres- electron cyclotron line features have been detected, suggesting 

^ ■ ence of neutron stars with strong magnetic fields. Observed spin that B ~ 10'^ - lO'^ G according to the formula for the elec- 

^ ; periods and their time derivatives in conjunction with the as- tron cyclotron energy, E^^^ HB/imeC) = ll.SSCB/lO'^ G) keV 

■ sumption of a magnetic dipole radiation give us magnetic field dOriandini & Fiume ||200l[). For many other X-ray pulsars with 

• !-( sti-ength as B oc (PP)^'^ where P and P are the spin period no detectable electron cyclotron line features, typical magnetic 

^ and its derivative, respectiv ely. For radio pulsars, of which more fields are B ~ lO'^ G if one assumes that spin u p due to accre- 

H ■ than 1 800 are known today ([Manchester et al.l2005h . the infeiTed tion o f matter is balanced by magnetic braking (jBildsten et alj 

... value of the magnetic field strength is in the range lO^'-lO'"^ G. |1997[). 

For a smaller population of older, millisecond pulsars, the typical From a theoretical point of view, these strongly magnetized 

magnetic field strength is B ~ lO^'-lO'' G. For anomalous X-ray objects deserve a detailed study, because their magnetic fields 

pulsars (AXPs) and soft gamma repeaters (SGRs), super strong sometimes exceed the quantum critical field strength, Bqed = 

magnetic fields of 10'^ - 10 '^ G are again inferred from the mea- m^c^/ieh) - 4.144 x 10'^ G, at which gyration radius of the 

sured values of P and P (Wood s & Thompson 11200 4'). Various electron pc/(eB) is shorter than the de Broglie wavelength h/p. 

observed properties of AXPs and SGRs Uke the giant flares from Above this limit, magnetic fields affect t he pro perties of atoms, 

the three SGRs and bursts are o ften explained in connection molecules, and condensed matters (lLaill200li) . propagation of 

with a supe r strong magnetic field (iThompson & Duncan II 19951 photons, radiative proce sses, equation of states, and thermal con- 

[T996ll200l rather than with rotation because their spin down lu- ductivity in crusts (see iHarding &Lail (I2OO6) and references 

minosities are much smaller than the observed luminosity. In ad- therein). The origin of such extremely large magnetic fields 

dition, temporary detections of spectral line s during SGR/AXP has been a big issue since their discovery. Specifically, there 

bursts have been reported in several sys tems (iGavriil et al. 120021: are two hypotheses for its origin; in one scenario, the mag- 

llbrahim et al.1l2003t iRea et al. II2003I) . If we assume that they netic fields are assumed to be generated in a rapidly rotating 
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prot o-neutron star formed after ste llar core collapse of a massive 
star (iThompson & Duncan 1ll993h and in the other, it is assumed 
to be descended from the main sequence stars, i.e., the strong 
magnetic field is ass umed to be a fossil of a strongly ma gnetized 
main sequence star dWickramasinghe & Ferrario 112005b . For ex- 
ploring the magnetar formation, there are a plenty of magneto- 
hydrodynamic simulations for supernova core collapse both 
in Newtonian gravity (Yamada & Sawai 2004; Kota ke et al. . 
2004; ObergauHngeretal. 2006; Scheideg geFet al. I 120081; 
Takiwaki et al. 2009; Burrows et al. 2007), and in ge neral 



relativity (Shibata et al. 2006; Cerda-Duran e t al. I l2007b . For 
these works, the simulations were performed in axisymmet- 
ric spacetime. In most of these simulations, it was found that 
toroidal magnetic fields are dominantly enhanced in the proto- 
neutron stars after the core bounce by the magnetic winding 
mechanism, and eventually, the proto-neutron star settles to a 
quasi-equilibrium state. However, it has been well known that 
purely toroidal fields in equilibria are often unstable due to 
the interchange, Tayler, and Parker instabilities dAchesonlll978t 
lGoossensllT980l: iParker IIT9551 IT9661 iTavler iflgTl . The stabil- 
ity analyses have suggested that non-axisymmetric modes would 
play an essential role for these instabilitie s and rapid rotation 
could suppress these instabilities (. Achesonll 1 9781) . R 

Motivated by these facts, we studied the axisymmetric insta- 
bility of neutron stars with toroidal magnetic fields in the pre- 
vious work CKiu chi et al. 2008). In that work, magnetized neu- 
tron stars in equilibria are prepared as initial conditions with 
varying its profile and strength and with changing the angu- 
lar velocity (Kiuchi & Yoshida 2008). Performing the general 
relativistic magneto-hydrodynamic (GRMHD) simulations, we 
found that slowly rotating neutron stars with the toroidal fields, 
whose profile is proportional to the power of cylindrical radius 
tJT as Z?(^) oc ?ir^*"' with k > 2, are unstable. The growth time 
scale is of order the Alfven time scale, and the type of the in- 
stability is the interchange instability in the absence of stellar 
rotation. Only for the case k = I, slowly rotating neutron stars 
are stable against the axisymmetric perturbation, and thus, we 
concluded that the configuration with k - 1 will be the attractor 
for the unstable neutron stars. We also found that rapid rotation, 
with which the rotational kinetic energy is much greater than the 
magnetic energy, suppresses the onset of the interchange insta- 
bility, and stabilizes the neutron stars. These results qualitatively 
and semi-quantitatively agree with the local linear analysis. 

In the magneto-rotational explosion scenarios, the magnetic 
energy is composed primarily of the toroidal field and is at most 
as large as the rotational kinetic energy. This indicates that the 
axisymmetric interchange instability would not play an impor- 
tant role in the proto-neutron star. However, it is still possi- 
ble that the neutron star becomes unstable against the Parker 
and Tayler instabilities which grow in a non-axisymmetric 
way. Stability of magnetars is also the important issue along 
this line (Braithwaite & Nordlund 2005; Braithwaite & Spruit 
12004 ). For the magnetar, the rotational effect is negligible. Thus, 
the toroidal field profile should be similar to that of A: = 1 
to avoid the onset of the axisymmetric interchange instability. 
However, such configuration may be still unstable if the Parker 
and/or Tayler instabilities are taken into account. 



' Besides the instabilities listed here, magneto-rotational instability 
(MRI) (Balbus & Hawlev 1991) could play an important role for en- 
hancing t he magnetic field strength and for modifying the magnetic field 
profile ( Obergaulinger et aTlllOOd IShibata et al. II2OO6V However, any 
accurate MHD simulation, in which the fastest growing mode of MRI 
is resolved, has not been performed yet. 



M otivated by these facts, we extend our previous 
work (iKiuchi et al.1 12008*). The main aim of this article is 
to explore the non-axisymmetric instabilities of neutron stars 
with toroidal magnetic fields. Following our previous work, 
we prepare equilibrium neutron stars with purely toroidal 
magnetic fields as initial conditions. This time, we perform 
three-dimensional GRMHD simulations with varying the field 
strength and/or rotation velocity. As mentioned above, the three- 
dimensional simulation is inevitable for exploring the Parker and 
Tayler instabilities. 

This paper is organized as follows. In Section |2] we briefly 
review the formulation and numerical methods employed in our 
numerical-relativity simulation. Set up of numerical simulation 
and initial models for our GRMHD simulations are described in 
Section |3] In Section |4l we present the results of a local linear 
perturbation analysis as a forecast of numerical-simulation re- 
sults. Section|5]is devoted to presenting the numerical results. A 
summary and discussion are given in Section |6] 

Throughout this paper, we adopt the geometrical units in 
which c = G = 1 with c and G being the speed of light and gravi- 
tational constant, respectively. Cartesian coordinates are denoted 
by X* = (x, y, z). The coordinates are oriented so that the rotation 
axis is along the z-direction. We define the coordinate radius 
r - -sjx^ + z}, cylindrical radius lu - + y^, and az- 
imuthal angle ip - tan"'(37/ji:). Coordinate time is denoted by f. 
Greek indices yU, v, ■ ■ ■ denote spacetime components, and small 
Latin indices /, j,--- denote spatial components. 



2. Formulation and Method 

The stability of magnetized neutron stars are studied by three di- 
mensional GRMHD simulation assuming that the ideal MHD 
condition holds. In this paper we focus on the Parker and/or 
Tayler instabilities against non-axisymmetric perturbations. The 
simulation is performed upgrading our jixisjinmetric GRMHD 
numerical code in IShibata & Sekiguchil (l2005i) : iKiuchi et aP 
(2008) to that for three dimensions. 

Formulation and numerical scheme for solving Einstein's 
equation are essentially the same as in Shibata & Nakamura I 
(1995). For solving Einstein's evolution equation, we use 
the original version of the B aumg arte- Shapiro - S hibata- 
Nakamu ra formula t ion dShibata & Nakamura I Il995l 
Baumgar te & Shapiro I Il999h : We evolve the inverse square 
of the conformal factor W = exp(-20) with (p - ln(y)/12, 
the trace part of the extrinsic curvature, K, the conformal 
three-metric, y,j = y^^^^jij, the tracefree extrinsic curvature, 
Aij = y^^^^(Kij — KjijlJ), and a three-auxiliary variable, 
Fi = d-'^'djfik- Here, y,y is the three-metric, Kij the extrinsic 
curvature, y = det(y,j), and K = Kgy'K Note that w e evolve 
W, not the conformal factor as in Kiuchi et al. ("20091), because 
our code is designed to si mulate bl ack hole spaceti mes with the 
moving puncture metho d ( Baker et al. II2006I: ICampanelli et al. I 



l2006l : lBrugmann et al. 1120081) . For the conditions of the lapse, a, 
and the shift vector, /?' we adopt a d ynamical gauge condition 
in the following forms (IShibata Il2003h . 



(5, -/3'5,)lnQ' = -2K, 
d,0 = Q.lSy'KFj + Atd,Fj), 



(1) 
(2) 



where At denotes the time step in the numerical simulations, and 
the second term in the right-hand side of (|2} is introduced for 
stabilizing the numerical computations. The finite-differencing 
schemes for solving Einstein's equation is essentially the same 
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as those in iKiuchi et al. I (l2009h . We use the fourth-order finite- 
differencing scheme in the spatial direction and a fourth-order 
Runge-Kutta scheme in the time integration, where the advection 
terms such as /S'djW are evaluated by a fourth-order non-centere d 
difference scheme, as proposed, e.g., in lBrugmann et al.1(l2008b . 

A conservative shock-capturing scheme is employed to 
integrate the GRMHD equatio ns. Specifically we use a 
high-resolution cent r al sch eme (iKurganov & Tadmor I 120001: 
iLucas-Serrano et ai~l |2004|) with the third-order piece-wise 
parabolic interpolation and with a steep min-mod limiter in 
w hich the limiter parameter b is set to be 2.5 (see appendix A 
of lShibatalJIOOl ). 

Magnetized neutron stars in equilibrium, employed as the 
initia l condition, are computed giving the polytropic equation of 
state (IKiuchi & Yoshida II2008I) . 



P = Kp^, 



(3) 



where P, p, k, and Y are the pressure, rest-mass density, poly- 
tropic constant, and adiabatic constant, respectively. In this 
work, we choose T - 2. Because k is arbitrarily chosen or else 
completely scaled out of the problem, we adopt the units of k - \ 
in the following (i.e., the polytropic units of q.=G=k - 1 are em- 
ployed). In the numerical simulation, we adopt the T-law equa- 
tion of state as 



p = (T-\)pe, 



(4) 



where s is the specific internal thermal energy. 

We monitor the total baryon rest mass Mt,, Arnowitt-Deser- 
Misner (ADM) mass M, internal energy f/int, rotational kinetic 
energy T^oi, total kinetic energy T^^n, and magnetic energy //mag, 
defined by 



(5) 



Mh - J" pwyfyd^x 

J e-^[phw^ -P+ -^{KijK'j -K^- Re''^''')] yfyd^ie) 



Uint - ^ pwsyfycf'x 

Tiot = ^ J phwu^Q.y/yd^x 

Tkin = ^ J phwuiv' ^Jyd^x 



(7) 
(8) 
(9) 
(10) 



where R is the Ricci scalar with respect to fjj, m'^ is the four 
velocity of the fluid, h is the specific enthalpy defined by 1 + 
s + Pip, w = au', v' = u'lu', Q. = vf, and b- = Vb^,. V is the 
magnetic field observed in the frame co-moving with the fluid 
element. Mo denotes the initial value of the ADM mass. Once 
each energy component is obtained, we define the gravitational 
potential energy by 

W = M - (Mfa + f/,„t + Tkin + //mag). (11) 

Following IKiuchi et alJ ll2008h . we define an averaged Alfven 
time scale as 



va = 



2H„ 



A|M/,+mmt+2//mag' 



(12) 



where we use the relation h - \ + Fe, which holds in the F-law 
equation of state. Note that the Alfven velocity in relativity is 



given by yjb^liAnph + b^). Then, we define the averaged Alfven 
time scale as 



R 

ta = —, 

Va 



(13) 



where R is the equatorial stellar radius. Because the magnetic 
field instability grows on an order of the Alfven time scale (cf. 
Section |4]i, ta is useful to judge whether or not the instability 
found in numerical simulation is associated with a magnetic field 
effect. 



3. Model and Numerical setup 

3.1. Initial condition 

Neutron stars with toroidal magnetic fields in equilibrium, em- 
ployed as initial conditions, are computed by the code described 
in Kiuchi & Yoshida (2008). Several key quantities character- 
izing these magnetized neutron stars are listed in Table |T] The 
instabihty associated with the presence of toroidal magnetic 
fie lds depends on t h e profile o f the magnetic field a s shown 
by AchesonI (Il978l) : iGoossens I (Il980l) : ISpruitI (Il999l) : iTavlerl 
(1973). We assume the toroidal magnetic field profile confined 
inside the neutron star to be given by 



(14) 



where k and Bo are constants which determine the field pro- 
file and strength, respectively. The regularity condition of mag- 
netic fields near the axis of zrr = requires k > \. Because 
of y^p^ oc up- for tir — > 0, the toroidal magnetic field is pro- 
portional to ziT^* ' near the axis. Previous works predict that 
the profile with ^ > 2 is unstable a gainst axisyrn metric pertur- 
bations (iAchesonlll97l iGoossens lll980t ISpruit 111 9991: iTavlerl 
197i and we confirm ed this prediction in the previous pa- 
per ( Kiuchi et al.ll2008h . On the other hand, the profile of ^ = 1 
is not unstable against axisymmetric perturbatio ns but may be 
unstab l e against non-axi s y mmetric ones (see also lLander et al. I 
(l2010h : lLander & Jones I (120101) ). Hence in this paper, we focus 
on the profile oik - 1 . 

Magnetic field strength, Bq, is chosen so as to satisfy 8 x 
10"^ < ///|W| < 5 X 10"-. These values imply the field strength 
of 10'*-10'' G for a typical neutron star of mass I.AMq, radius 
10 km, and \ W\ ~ 6 x 10^^ erg. These magnetic field strengths 
are extremely large even for models of magnetar and might be 
less realistic. We here give such strong magnetic fields simply to 
save the computational costs; note that the growth time scales of 
the instabilities by the presence of the toroidal magnetic field are 
of order Alfven time scale, which is still much longer than the 
dynamical time scale of the system in the present set up. Thus, a 
scaling relation should hold for a weaker magnetic field strength. 
We may apply the scaling relation to derive a generic physical 
essence from the results obtained in the present set up. Namely, 
if the magnetic field strength becomes half, the growth time 
scale of the instabilities becomes approximately twice longer, 
although qualitative properties of the evolution of the unstable 
neutron star are essentially the same. 

The initial conditions for the non-rotating model is specified 
if the central density pc (in the polytropic units) is determined. 
We choose it to be p^ ~ 0.22. With this choice, the neutron star 
has a realistic compactness; e.g, if we assume M » 1.35M0, 
circumferential radius is R ^ 11km (cf. Table [T]i. Note that the 
maximum rest mass and gravitational mass of the spherical neu- 
tron star for F = 2 are, respectively, 0.1799 and 0.1637, and the 
corresponding central density is » 0.318 in our unit. 
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Simulations are also performed for rotating neutron star 
models. In this paper, we focus only on rigidly rotating neu- 
tron stars with a moderate compactness; pc is again chosen to 
be 0.22. For the neutron stars with the F = 2 polytro pe, the max- 
imum values of r,-ot/| W| is ~ 0.09 for Mq //^cir ~ 0. 1 dCook et al. I 
|1994|) . Local linear stability analysis predicts that rapid rotations 
will suppress the growth of some of the instabilities associated 
with the presence of toroidal magnetic fields (Acheson (1978), 
see also Section 4.2). To study whether this would be indeed the 
case, we prepare rapidly rotating models with Trot/IWI ^ 0.08, 
and vary the magnetic field strength. 



3.2. Grid settings 

The simulations are performed on the cell-centered Cartesian, 
(x, y, z), grid. Equatorial plane symmetry with respect to z = 
plane is assumed. The computational domain of -L < x < L, 
—L < y < L, and < z < L is covered by the grid 
size (2N, 2N,N) for (x,y,z), where L and are constants. 
Following iKiuchi et alJ (l2008h . we adopt a non-uniform grid as 
follows; an inner domain is covered with an uniform grid of 
spacing Ax and with the grid size, {2No,2N(),Nq). Outside this 
inner domain, the grid spacing is increased according to the re- 
lation, ^tanh[(/ - N{))/Ai]Ax, where / denotes the i-th grid point 
in each positive direction, and A^o, A/, and ^ are constants. Then, 
the location of i-th grid, x''(i) (i > 0), for each direction is 



, , (/+l/2)Ax '-''h5) 
^ \ a + l/2)Ax + ^MAx log[cosh{(/ - M))/A/)] / > m ' 



i < Na 



and x*(-/ - 1) = -x''(i), where / = 0, 1, ■ ■ ■ - 1 for x*^ - x,y, 
and z. The chosen parameters of the grid structure for each sim- 
ulation are listed in Table|2] To check the validity of our numeri- 
cal results, the simulations for models N22H5 and R22H2T8 are 
performed with three different grid resolutions. With the high- 
est (lowest) resolution, the coordinate equatorial radii of neutron 
stars, Lns, are covered by 130 (80) grid points while the outer 
boundary location is chosen to be ^ 8Lns in all the simulations. 
We confirmed the modest resolution in which Lns is covered 
by 100 grid points is high enough to draw a resolution-invariant 
conclusion. 



4. Local linear perturbation analysis 

Before showing the results of nonlinear numerical simulation, 
it is useful to remind the results of linear- stabihty analysis . By 
using the local linear perturbation analysis (lAcheson Ill978h . the 
dispersion relations for the non-rotating and rotating models are 
derived, and the (local) stability is determined. In the following, 
the method of the analysis in the Newtonian framework is briefly 
reviewed. 

In the local analysis, a perturbation quantity 6Q is assumed 
to be given by 



= Qo exp[/(/tir -I- m(p + nz- erf)]. 



(16) 



where Qq is a constant, cr is the oscillation frequency, and 
(/, m, n) is the wave number vector. To obtain the local dispersion 
relations, we assume that / » 1, n » 1, and m - 0(1), which im- 
plies that perturbations we consider ha ve very short wav elength 
on the meridional plane (for details, cf . lAcheson I (1 1 9781) ) . 



4.1. Non-rotating case 

The dispersion relation for non-rotating models is given, irre- 
spective of the density profile, by 



.4 



2 2 

m cr. 



P 



+ 1 + ,„V^ + m-'crlm' \^-—\dh \n(B^^^m) = (17) 



■UT 



where va = ^J4JTp, cta = va/tu, G ^ g,^ - (lln)g, with g,^ 
and g~ being the gravitational acceleration, 5/, = 5nj - {lln)d,, 
and Cj is the sound speed. For the axisymmetric perturbation, 
i.e., m - Q perturbations, the dispersion relation is reduced to a 
quadratic form of cr as 



.4 

n- 



1 + -4- 



= 0. (18) 



For the density profile which is necessary to solve Equation 
(fTTl l for a specific model, we prepare the Newtonian spherical 
polytrope with F = 2 as 



P =Po- 



s in(7rr/r^) 
nr/r. 



(19) 



where po and are the central density and stellar radius. The 
spherical density profile is justified by the assumption of a weak 
magnetic field, i.e., the magnetic field is too weak to deform the 
star. In the Newtonian limit. Equation (fT4l l becomes 



(20) 



Because B^^)/pm =const. in our models (see Equations ( ST% 
and (l20li). we have cr - Q for m = 0, i.e., the magnetic 
field is marginally stable against the axisymmetric perturbation. 
Physically, the magnetic field instability caused by the axisym- 
metric perturbation is classified as the interchange instability. 
Our models are marginally stable against such instability and 
we i ndeed found in the previous study that they were not unsta- 
ble (IKiuchi et alJl2008h . However, the marginally stable profile 
is not always kept in a stationary state in the presence of a per- 
turbation. We return to this point in Section 15721 

For non-axisymmetric perturbatio ns, we forecast that the 
Par ker instability (|Parker1ll955l 1 19661) and the Tay ler instabil- 



ity ( Tavler II 19731) set in. Following Acheson! (1 19781) . we define 
the critical radius as 



2c 



(21) 



Outside this critical radius, the magnetic buoyant force ex- 
ceeds the magnetic hoop stress, whereas inside it, the magnetic 
stress is dominant . The Tayler instability could set in near the 
axis of T3T = (iSt?ruitlll99 9^. which is inside the critical ra- 
dius. On the other hand, the Parker instability could play a role 
in places where the magnetic buoyancy force surpasses the mag- 
netic tension. Thus, we classify the instability which emerges 
inside and outside the critical radius as the Tayler and Parker 
instabilities, respectively. Analyzing the dispersion relation near 
the axis of tir = 0, Acheson ( 1978) and Spruit ( 1999) concluded 
that the dominant modes of the Tayler instability are m - 1 and 
l/n X Q modes, which are associated with the horizontal motion 
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Table 1. List of characteristic quantities for neutron stars with toroidal magnetic fields. 



Model 




C/m,/|W| 




T,J\W\ 


Mo 


Mb 




ta/Mo 


N22H5 
N22H1 


2.2E-1 (1.24) 
2.2E-1 (1.24) 


5.12E-1 
5.32E-1 


5.01E-2 
l.OlE-2 






1.59E-1 (1.67) 
1.60E-1 (1.68) 


1.73E-1 (1.81) 
1.75E-1 (1.96) 


1.78E-1 
1.87E-1 


37.4 (0.31) 
76.7 (0.63) 


R22H2T8 
R22H08T8 


2.2E-1 (1.24) 
2.2E-1 (1.24) 


4.7 IE- 1 
4.74E-1 


1.50E-2 
7.86E-3 


7.97E-2 
8.10E-2 


1.84E-1 (1.93) 
1.85E-1 (1.94) 


2.01E-1 (1.63) 
2.02E-1 (1.73) 


1.55E-1 
1.65E-1 


63.5 (0.63) 
97.3 (0.93) 



Notes. Central density (p^), ratio of the internal energy to the gravitational potential energy W (t/i„t/|W|), ratio 
of the magnetic energy to W (//mas/IM^I), ratio of the rotational kinetic energy to W (T^ot/\W\), ADM mass (Mo), 
baryon rest mass (Mt,), and compactness (Mo/R^-iy) with iJ^ii being equatorial circumferential radius, is an 
averaged Alfven time in units of Mo. Model name NXHY denotes non-rotating models with "X" and "Y" being 
the values of lOOpr and the values of 100//„,a8/|W'L respectively. RXHYTZ denotes rotating models with "Z" 
being the values of 100r,ot/| W|, where meanings of "X" and "Y" are the same as the non-rotating models. The 
values shown in the brackets denote those in physical units, where the densities are normalized by lO'^g/cm'', 
the masses by M©, and the Alfven time are given in units of millisecond. Note that we set as 1.6 x 10^ in cgs 
units. 



Table 2. Parameters for the grid structure employed in the numerical simulation. 



Model 


N 


A'o 


M 




Lns/Ax 


N22H5-1 


161 


120 


30 


22 


80 


N22H5 


200 


150 


30 


20 


100 


N22H5-h 


253 


195 


30 


20 


130 


N22H1 


200 


150 


30 


20 


100 


R22H2T8-1 


161 


120 


30 


22 


80 


R22H2T8 


200 


150 


30 


20 


100 


R22H2T8-h 


253 


195 


30 


21 


130 


R22H08T8 


226 
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Notes. The grid number for covering one positive direction (N), that for the inner uniform grid zone (A'o), the 
parameters for nonuniform-grid domain (Ai, f ), and the approximate grid number for covering the coordinate 
radius of neutron star (Lns), respectively. 



near the axis of zir = 0. However, to clarify the most relevant 
unstable modes, we have to determine the growth time scales of 
the instability for all the places inside the star. Hence, we cal- 
culate the maximum growth rates of the instability solving the 
dispersion relation (fTTl i with varying l/n and m. 



4.2. Rotating-case 

For the rotating models, the dispersion relation is written as 
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(22) 



Figure [T] displays the contours of the growth rate for the 
model with po - 0.16 and Bq - 0.2, which give H^^g/\W\ - 
2.5% and Mi, = 0.2 in the polytropic units. These parameters 
are chosen so as to mimic the initial models shown in Table [1] 
The white colored region and the dotted curve denote the stable 
region and critical radius, respectively. We note that the most un- 
stable mode at each point is the l/n - mode. This figure shows 
that the fastest growing mode of the instability is located near 
the stellar surface and is determined by the Parker instability, be- 
cause the location is outside the critical radius. The right panel 
of Figure [1] shows that not the m = 1 mode but a high-order m 
mode is relevant for this fastest growing mode. The local linear 
perturbation analysis predicts that the Parker instability primar- 
ily emerges near the stellar surface in our non-rotating models. 



where to — cr - otQ w ith Q being the angular velocity. 
Following lAchesonl(ll978h . we focus on a low-frequency and 
non-axisymmetric mode, for which the growth time scale is 
much longer than the Alfven time scale as 

\Lof- <K mV^. (23) 

We also adopt a weak magnetic field approximation in which 
magnetic energy is assumed to be everywhere much smaller than 
the rotational kinetic energy as 

vl <K n^m^. (24) 
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With these approximations, the bi-quadratic equation 
duced to a quadratic equation and its solution is 



IS re- 



'G_ 



CJ^\^ 



mvA 



ml 2QcT 2Q.TIT 



^ - — I v\dh In 



\pmj 



2 2 



1/2 



(25) 



This dispersion relation shows that the instability sets in if either 
of the following criteria is satisfied, 



cj \pmj 



2 2 
> m o- . 



(26) 
(27) 



The inequalities ( l26b and (l27b do not hold for rigidly rotating 
stars with the magnetic field profile (l20t . However, the angular 
velocity profile of any star never remain constant in a strict sense, 
because the stars in nature usually oscillate around their equilib- 
ria which makes a gradient in the angular velocity profile. In the 
presence of negative angular velocity gradient, the instability cri- 
terion (|26] | may be satisfied because the first term in the left-hand 
side of Equation (|26] | can be the most dominant one due to the 
weakness of the magnetic field as imposed in Equation (l24l i: for 
a region in which the Alfven velocity is much smaller than the 
rotational velocity, a small perturbation in the angular velocity is 
sufficient for satisfying the inequality (l26T l. We will discuss this 
issue in Sectionl572lin more detail. 



5. Result 

5.1. Non-rotating case 

We study the stability for two non-rotating models listed in 
Table [T] The simulations are performed for a sufficiently long 
time, more than ten times of the averaged Alfven time scale or 
several thousands of time in units of Mq (cf. Table [T]). This is 
necessary to clarify whether or not any MHD instability, which 
grows approximately on an Alfven time scale, sets in and to de- 
termine the final fate after the onset of the instabilities. Figure |2] 
plots the evolution of the central rest-mass density and the min- 
imum value of the lapse function, which characterize the com- 
pactness of the neutron star. For the non-rotating model N22H5, 
it is observed that the central density (the minimum value of the 
lapse function) increases (decreases) for t < 200Mo and subse- 
quently decreases (increases) for 200 Mq < t < 600M() and then 
stably oscillates for t > 600M(). This suggests that the star con- 
tracts first, then expands slightly, and finally settles to another 
quasi-equilibrium state with oscillations. Note that the final cen- 
tral density is not beyond the marginally stable point, X 0.318. 
For model N22H1, the qualitative features of the evolution are 
essentially the same as those of N22H5, but the time scale is 
different from that of N22H5 because of the difference in the 
magnetic field strength. The behavior described above is well 
explained by the variation of magnetic fields during the evolu- 
tion as argued below. 

Figures [3] and |4] plot the evolution of the rest-mass density 
and magnetic energy density on the equatorial plane and in one 
of meridian (x-z) planes, respectively, for model N22H5. The 
panels (a)-(c) in these figures show that the magnetic field near 
the stellar surface is disturbed by the Parker instability, and leaks 
out of the stellar surface. Because the plasma beta is small and 



thus the matter inertia is small near the stellar surface (as shown 
below), the matter is dragged by the magnetic force and conse- 
quently the stellar surface is distorted. Here, the plasma beta is 
the ratio of the fluid pressure to the magnetic pressure. 



= ^ A 1_ 

^p,.™a = oc ^^^^ oc 



(28) 



where we have used Equations Q and ( fT4l i assuming F = 2 
and ^ = 1 . The minimum value of the plasma beta is initially 
^ 2 at the stellar surface, and after the onset of the Parker in- 
stability, the leak-out magnetic field loop produces even lower 
beta plasma near the stellar surface. As a result, a weak wind ex- 
panding outward is driven. On the other hand, the ingoing mag- 
netic field loop enhances a turbulent motion in the neutron star 
During the transition from the state shown in panel (c) to (d) in 
Figures [3] and |4] the initial magnetic field profile is completely 
destroyed and turbulent magnetic field is produced. During the 
development of the turbulence, the Tayler instability does not 
appear to play an important role. However, the region near the 
axis of t(7 = is not stable against this instability and thus no 
mechanism seems to help stabilizing there. 

The toroidal magnetic fields initially prepared behave like a 
rubber belt, which fastens the "waist" of neutron stars. The dis- 
appearance of the coherent toroidal magnetic fields, therefore, 
results in the expansion of the star as shown in the panel (d) 
of Figures [3] and |4] After the magnetic field becomes turbulent, 
the star stably oscillates around the hypothetical quasi-stationary 
state. Although the density profile relaxes to a quasi-stationary 
state, the turbulent motion is maintained. We observe qualita- 
tively the same features for model N22H1, but the growth time 
scale of the instability is longer than that of N22H5 as mentioned 
before. 

It is interesting to compare the simulation result with the lin- 
ear analysis. We find the turbulent field develops in the region 
which is stable against the perturbation (see the stable region 
in Figure[T]and Figure|4](d-f)). During the linear growth phase, 
i.e., Figure|4](a-c), this region is likely to be stable. However, be- 
cause the instability destroys the coherent initial magnetic field 
profile, i.e., the magnetic field is no longer pure toroidal, the sta- 
ble region in Figure [T] is no longer stable after the onset of the 
instability. We also point out that the linear analysis in Figure [1] 
indicates that the higher m mode is unstable near the stellar sur- 
face and we indeed find such a behavior (see Figure |3](b)). For 
making the mode growth clear, we perform the mode analysis as 
follows. Because we are interested in the magnetic field behav- 
ior near the surface, we put a ring on the equator whose radius is 
nearly equal to the stellar radius. Then, the Fourier components 
are defined by 



C,. 



-i ' 



(29) 



where r^ng is the ring radius which is chosen as AMq for model 
N22H5. Figure|5](a) plots the evolution of |C,„| with 1 < m < 10. 
It is found that all the mode (except for m - 4 and 8) start grow- 
ing at r a; 150Mo, that shows that the Parker instabilities are in 
operation. For m - A and 8, \Cm\ starts growing at T a; 50Mo. 
This is the artifact due to our choice of the Cartesian coordinate 
grid. 

The nonlinear growth of the Parker instability is well re- 
flected in the energy components as shown in Figure |5] Here, 
the internal energy ^ is separated into the adiabatic part t/ad 
and the heating part C/heat, and they are defined by replacing 
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E in Equation d?) to Sad = V(r - 1) and to e - Sad, re- 
spectively. Figure |5] shows that the ADM mass is approximately 
conserved, that implies that gravitational waves are not substan- 
tially emitted during the evolution of the neutron stars. The mag- 
netic energy gradually decreases and the kinetic energy increases 
during t < 500Mo for model H22M5. This illustrates that the 
meridional circulation is excited due to the instability: Figures[3] 
and ID specifically show that this meridional circular motion is 
induced by the displacement of the toroidal magnetic field line. 
In other words, the magnetic field helps increasing the kinetic 
energy of the fluid element by liberating the gravitational poten- 
tial energy. During this procedure, the kinetic energy eventually 
reaches about ten percents of the magnetic energy at f =s 200Mo. 

For t > 500M(), the decrease rate of the magnetic energy 
becomes high and simultaneously the thermal energy (t/heat) 
quickly increases. This indicates that the shock heating oc- 
curs because of the turbulent motion induced by the distorted 
magnetic fields. The kinetic energy of the circulation gradu- 
ally decreases, and the magnetic and thermal energies settle 
approximately to constant values at a late stage of the evolu- 
tion. The adiabatic internal energy also settles to a value lower 
than the initial one at f a; SOOMq. This implies that the den- 
sity distribution changes, as shown in Figures |3] and |4] As we 
pointed out above, however, the magnetic field configuration 
never reach es to any equilibrium state in contrast to the axisym- 
metric case dKiuchiet al.ll2008h . 

All the features found in the evolution of the energy com- 
ponents for model N22H1 are essentially the same as for 
model N22H5 except for the relevant time scale, as found from 
Figure s|2] and |5] The variation time scales of the various compo- 
nents of the energy for model N22H5 are systematically shorter 
than those for model N22H1. This is natural because the insta- 
bility develops approximately in the Alfven time scale. 

As found from Figure |4] the meridian circulation is excited 
by the magnetic field instability. Thus, the growth time scale of 
the kinetic energy in the rising-up phase (f < IQQMq) should be 
of order of the Alfven time scale. To check if this is indeed the 
case, we evaluate the growth time of the kinetic energy. For this 
evaluation, the data of f/M,, e [50, 100] and of f/Mo e [50, 200] 
are used for models N22H5 and N22H1, and are fitted assuming 
the function form of oc e'/'"' where is the growth time. We find 
that Tk ~ 14Mo and ^ 22Mq for models N22H5 and N22H1, 
respectively. Because these values are substantially smaller than 
the averaged Alfven time scale given in TablelT] ta is less appro- 
priate for characterizing the growth time scale of the magnetic 
field instability. Rather, we find it appropriate to employ Alfven 
time scale in the vicinity of the stellar surface because both the 
linear analysis and our simulations indicate that the instability 
primarily grows there. The Alfven time scale ta estimated at 
r/r, = 0.95 is 60Mo for model N22H5 and 120Mo for model 
N22H1. Furthermore, the linear analysis suggests that the mode 
of m ~ (9(10) is most unstable and we find such a feature in 
our simulations as discussed above. Because the growth rate is 
proportional to m as we see in Equation (ITTT i. the growth time 
scale should be given by TA/in, ~ 6Mo for model N22H5 and 
12Mo for model N22H1, which agrees with within a factor 
of two. In any case, the growth time scale of the instability is 
approximately proportional to the Alfven time scale normalized 
by m. Therefore, we can conclude that the primary instability is 
the Parker instability as expected in the linear analysis. 



5.2. Rotating case 

We study the stability for two rigidly rotating models listed in 
Table [T] Again, long-term simulations are performed as in the 
non-rotating models. In Figure |2] the evolution of the central 
density and the minimum value of the lapse function for models 
R22H2T8 and R22H08T8 are plotted together. We note again 
that for both models, the angular velocity is approximately as 
large as the Kepler limit at the equatorial stellar surface. The 
central density for model R22H2T8 (R22H08T8) remains ap- 
proximately constant until t » 800M()(1000Mo), and then, grad- 
ually increases for t < 2000Mo(2400M()). Eventually, it settles 
to a new value for t > 2000Mo(2400Mo). The reason that the fi- 
nal central density is larger than the initial one will be described 
below. 

Although the magnetic field strength and central density for 
the rotating model R22H08T8 are comparable with those for the 
non-rotating model N22H1, the evolution process in the cen- 
tral region is significantly different. For model R22H08T8, the 
central density remains approximately constant for a time much 
longer that for model N22H1. This is due to the fact that the ro- 
tation stabilizes the Tayler instability which may set in for the 
non-rotating models, as expected in the local linear perturbation 
analyses (see Section I4.2l i: The Tayler instability is a primary 
magnetic field instability associated with the toroidal field near 
the axis of tit = for the non-rotating models, but this is not the 
case for the rapidly rotating models. 

However, the instability is not suppressed by the presence 
of rapid rotation for all the places inside the neutron star, as 
Figure |2] illustrates that the central density varies for a longer- 
term evolution with t ^ 2000Mo. This result is totally different 
from that in the axisymmetric case (iKiuchi et al.l2008l) . in which 
the rapid rotation suppresses the onset of the interchange insta- 
bility. Note that one of the models studied in Kiuchi et al. (2008|) 
has similar model parameters to those of R22H08T8, which is 
stable for the axisymmetric perturbations. This implies that a 
non-axisymmetric instability is excited for the rapidly rotating 
models this time. To clarify the relevant instability, we generate 
Figures|6]and|7] in which the rest-mass density and the magnetic 
energy density for model R22H2T8 are plotted on the equato- 
rial and meridian planes, respectively. Until t ^ 200Mo, we ob- 
serve that the coherent magnetic field and density profiles re- 
main. However, at f ^ 400Mo, the magnetic field profile near 
the stellar surface is deformed in the same manner as found 
for the non-rotating models. During the subsequent evolution, 
the coherent magnetic field structure is totally destroyed and a 
turbulent motion is excited. The surface expands because the 
plasma beta near the surface is below unity due to the leak-out of 
the magnetic field. Because the coherent toroidal magnetic field, 
which fastens the waist of the neutron star, disappears, the radius 
of the neutron star increases. All these features are essentially the 
same as for the non-rotating models. 

The evolution of the magnetic field near the rotation axis 
in the rotating models is slightly different from that in the non- 
rotating models. Figures |8] (a) and (b) plot the snapshots of the 
magnetic pressure defined by l%n along the x and y-axes on 
the equator for models R22H2T8 and N22H5, respectively. In 
the rotating model R22H2T8, the profile near the center does 
not drastically change. However, the magnetic pressure increases 
gradually and systematically, and by the increased magnetic 
pressure, the matter near the axis is pinched. This increase pro- 
ceeds in the axisymmetric way because the profiles along y- 
axis is approximately identical to those along x-axis in Figure |8] 
(a). This might seem to be incompatible with our previous re- 
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suits (iKiuchi et alJ200^ . To clarify this point, we reconsider the 
criterion for the axisymmetric perturbation : 
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The axisymmetric instability ignites if any of these inequalities 
holds. Equations ( |3T] ) and ( l32b implies the magnetic field profile 
given by Equation (l20l i is marginally stable. The magnetic field 
line strength evolves as 



(5, + e ^ 'v(a)di) = 5,Q, 

\mp p 



(33) 



where e''"' is a tetrad basis. In axisymmetric simulations, the 
field line strength is completely preserved because the poloidal 
magnetic field, i.e. B^^-, and is never generated if the ini- 
tial magnetic field is purely toroidal. Then, the inequalities (ISTT i 
and ( l32l l are never satisfied during the evolution of the mag- 
netized star On the other hand, the field line strength deviates 
from its initial value in three-dimensional simulations because 
the poloidal field is generated by a perturbation induced by a 
truncation error in general (in addition to the change of the an- 
gular velocity profile as discussed below). As the result, the in- 
equalities (ISTT l and (l32l l may hold during the evolution and the 
axisymmetric instability could set in. This is the reason why 
the magnetic pressure increases and angular velocity is distorted 
near the center as seen in Figures |8](a) and (c). Note that the 
things are the same in model R22T08T8. It is interesting to note 
that the stellar surface in the rotating model expands whereas the 
central part contracts due to the redistribution of the magnetic 
field profile. Contrary to the rotating model, the coherent profile 
of the magnetic fields near the stellar center in the non-rotating 
model N22H5 disappears after the onset of the instability, which 
results in the systematic expansion of the star as mentioned 
in the previous subsection. The non-rotating models are also 
marginally stable against the interchange instability as shown in 
Equation ( ISTT i. However, these models are unstable against the 
non-axisymmetric mode (see Figure [T](b)) and this mode over- 
comes the the interchange mode. It is also interesting to note 
that the non-axisymmetric instability develops in an earlier time 
than that expected from the central density's evolution, e.g., at 
t * 800Mo(1000Mo) for R22H2T8 (R22T08T8) (see Figure |2] 
(a)). However, the instability develops before t ~ 6OOM0 for 
these rotating models as shown in Figure|9] This also reflects the 
fact that the instability, which affects the global structure of the 
neutron star, sets in both near the stellar surface and in the cen- 
tral region. Namely, the Parker (interchange) instability plays an 
important role in the outer (central) region for the rotating mod- 
els. 

The local linear perturbation analysis predicts that our rigidly 
rotating models are stable against the non-axisymmetric pertur- 
bation as described in Section 14.21 Thus, our numerical result 
does not seem to agree with that in the linear analysis. However, 
this is not the case, because our rotating model is close to a 
marginally stable state and is destabilized by a slight nonlinear 
perturbation to the rotational velocity. To explain this fact, we 
plot the snapshots of the angular velocity in Figure |8] in which 
the angular velocity profiles along the x-axis on the equator are 



plotted for models R22H2T8 and R22H08T8. This shows that 
the angular velocity deviates slightly from the constant profile 
and the negative gradients, in particular, near the stellar surface 
are developed during the evolution. This time variation is due to 
the oscillation of the neutron star which is initially triggered by a 
perturbation of numerical origin. Note that the initial conditions 
we gave are in equilibria, but a small numerical error associated 
with the finite grid resolution induces a perturbation and then 
the neutron stars start oscillating around their equilibrium states. 
Although the perturbation is induced by a numerical error in this 
case, it is quite natural to expect that any star in nature oscillates 
and thus the precisely rigid rotation is not realized. Once the an- 
gular velocity profile has the negative gradient, the instability 
criterion (l26T l could be satisfied because the first term in the left- 
hand side has a substantially large positive value in the rapidly 
rotating models with weak magnetic fields, i.e., the criterion is 
satisfied even for a small angular velocity gradient. 

The growth process of the instability is well reflected in sev- 
eral energy components. Figure |9] plots their evolution. This 
shows that the ADM mass, adiabatic internal energy, and ro- 
tational kinetic energy remain approximately constant during 
the evolution. The decrease in the magnetic energy is promi- 
nent at t K I5OOM0 for model R22H2T8 and t * 2000Mo 
for R22H08T8. Note that the magnetic energy increases only 
by ~ 2% during the development of the interchange instability, 
e.g, for t < 4OOM0 for model R22H2T8. The kinetic energy, in 
which contribution of the rotational kinetic energy is excluded, 
increases with time. This shows that the meridian motion is de- 
veloped. We find again that the rapid rotation is not enough to 
stabilize the instability associated with the presence of toroidal 
magnetic fields and the magnetic field never reaches an equilib- 
rium state as in the non-rotating models. 

It is well known that negative gradient of angular ve- 
locity in the presence of magnetic field leads to the 
MRJ (Balbus & Hawlev 1991). As discussed above, the nega- 
tive angular velocity appears in the vicinity of the stellar sur- 
face. Hence, we hypothetically estimate the MRJ growth rate as 
(Tmri ~ D.m\d^lnQ.\. We fit the angular velocity profile as a 
function of the radius in Figure |8] (c-d) and obtain the gradient 
c?cj In Q.. Putting the value of the gradient, angular velocity, and 
the stellar radius into the equation, we estimate the growth rate as 
o-MRiMo ~ 0.09 (0.1) for model R22H2T8 (R22H08T8). On the 
other hand, we infer the instability growth rate from the increase 
of the kinetic energy shown in Figure |9] with the same strategy 
for the non-rotating case. The resulting growth rate crMo is 0.008 
(0.007) for model R22H2T8 (R22H08T8), where we use the data 
of t/Mo e [0, 600] in both the models. Therefore, the growth 
rates do not match. However, Equations (fTTI i and (l25l l tell us that 
the growth rates of the interchange and Parker instability depend 
on the magnetic field strength. Hence, we conclude that the in- 
stability we find in the rotating model is not the MRI. 

Before closing this section, we show that the convergence is 
achieved in our numerical results with an accuracy high enough 
to draw a quantitatively reliable conclusion. Figures [TO] and [TTI 
plot the evolution of (a) the central density, (b) the minimum 
value of the lapse function, (c) the ADM mass, (d) the inter- 
nal energy, (e) the magnetic energy, (f) the kinetic energy, and 
(g) the Hamiltonian constraint violation for models N22H5 and 
R22H2T8. The Hamiltonian constraint violation is defined by 



ERROR =— f pw^/j\V\d\ 
Mb J 



(34) 
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with 

V = — — o_v_j J /_ ^22^ 

\Ai(r\ + if^^l + |27rpH<A'l + ^ {AijA'J + f/T^) 

where A is the Laplacian associated with y,y, tfr is the confor- 
mal factor, and pH - (ph + b^/47T)w'^ -(P + b'^l^n) - {aWflAn. 
In all the panels, we plot the results of three different grid res- 
olutions. In Figure [To] (a), (b), (d), and (e), we find a transition 
time dX t ~ 500Mo depends slightly on the grid resolution, but 
on the whole, the results appear to be convergent. The ADM 
mass, which is approximately conserved in the present situation, 
is conserved within 99.8 percent accuracy (see Figure [TO](c)) and 
the error in the Hamiltonian constraint violation is in an accept- 
able level (less than ~ 1% error; see Figure [Tn](g)). In the panel 
(e), we observe that the convergence is lost for t > IOOOMq. This 
is likely to be caused by the turbulent field developing both out- 
side and inside the star as shown in Figures |3] (f) and|4](f). To 
handle such a region and/or turbulent field, we would need more 
sophisticated numerical scheme. 

From these figures, we conclude that all the qualitative fea- 
tures of the magnetic field instability found in this work are inde- 
pendent of the grid resolution. In the rotating model R22H2T8, 
we find essentially the same features in the convergence study. 

6. Summary & Discussion 

6.1. Summary 

We explored the non-axisymmetric instability of neutron stars 
with purely toroidal magnetic fields. Preparing the non-rotating 
and rotating neutron stars in equilibrium as the initial conditions, 
the three-dimensional GRMHD simulations were performed. 
For the non-rotating models, the local linear perturbation analy- 
sis predicts that the Parker instability would be the primary in- 
stability and we confirmed this. Due to the Parker instability, a 
turbulent state is developed and the initially coherent magnetic 
field profile is totally varied. The magnetic field profile never 
reaches an equilibrium state. This fac t is in sharp contrast with 
that in the axisymmetric instability of iKiuchi et alJ (1200^ . The 
growth time scale of the Parker instability depends on the mag- 
netic field strength, i.e., the Alfven time scale, and this result also 
agrees with the local linear perturbation analysis. The present re- 
sult strongly suggests that three-dimensional treatment is crucial 
to clarify the instability of a neutron star with toroidal magnetic 
fields. In other words, any a priori assumption of the spacetime 
symmetry (e.g., axisymmetric symmetry) could prevent from de- 
riving the correct conclusion. 

We also explored the instability of rigidly and rapidly rotat- 
ing neutron stars. The linear analyses have suggested that rapid 
rotation could play a role as a stabilizing agent. We confirm that 
the rapid rotation stabilizes the Tayler instability, which may oc- 
cur near the axis of cr = in the non-rotating case. However, 
the interchange instability could play a minor role because the 
neutron stars are marginally stable against it. We note that the 
interchange mode never develops in the axisymmetric simula- 
tions due to the symmetry imposed. We also find that the Parker 
instability which is relevant near the stellar surface may not be 
stabilized by the rapid rotation. The reason is that by a pertur- 
bative oscillation, neutron stars may have a region in which the 
gradient in the angular velocity profile is negative (dQ./dm < 0). 
This negative gradient can induce the Parker instability in the 
case that the neutron star has rapid rotation and weak magnetic 



fields. As in the non-rotating model, a turbulent state is subse- 
quently developed in the outer region of the neutron star. This 
result also gives us a message that the three-dimensional simu- 
lation is essential for investigating a magnetic field instability. 

6.2. Discussion 

As mentioned in Introduction, a large number of core-collapse 
supernova simulations based on the magneto-rotational mecha- 
nism have shown that toroidal magnetic fields are dominantly 
amplified in the proto-neutron stars via the winding mechanism. 
In the assumption of axial symmetry, these fields may be in 
quasi-equilibrium after the saturation is reached. However, this 
could disagree with the result in the linear perturbation analy- 
sis, i.e., the neutron stars with purely toroidal magnetic fields 
are often unstable. The present study indeed suggests that such 
neutron stars are unstable and thus the assumptions of axial sym- 
metry and rapid rotation, which are imposed for most of the 
magneto-hydrodynamical supernova simulations, would be in- 
appropriate. (Note that in axial symmetry, the rapid rotation sta- 
bi lizes neutron s tars w ith toroidal magnetic fields, as illustrated 
in iKiuchi et al.l ( l2008l) .) In a non-axisymmetric simulation, we 
may find that a proto-neutron star with strongly toroidal mag- 
netic fields is unstable and a turbulent motion inside it is ex- 
cited. Magneto-hydrodynamic simulations have to be performed 
in fully three spatial dimensions. 

Stability of more generic magnetic field configurations, 
i.e., mixed poloidal-toroidal fields, is quite important be- 
cause such fields woul d be i n general realized. Recently , 
'Braithwaite & NordlundJ (l2005h : iBraithwaite & Spruit I i2004 : 
Duez et al. ( 2010) studied the stability of Newtonian stars with 
such fields, and reported that the mixed field may play an im- 
portant role in stabilization; they showed that the equilibrium 
profile is maintained over several Alfven time scales. We plan to 
stu dy this issue w i th the weak magnetic field solution obtained 
by iloka & Sasaki I (l2004l) in the fully general relativistic frame- 
work. 

Recently, 'Lander e t aP (1201 Oh and iLander & Jones I (1201 Oh 
studied the instability associated with the presence of toroidal 
magnetic fields by solving the linearized Newtonian MHD equa- 
tions with their time-domain code. They showed that the Tayler 
instability characterized by the azimuthal mode number m = 1 
primarily occurs near the axis of w = and that no pronounced 
Parker instability sets i n near the surface o f the star . Similar re- 
sults were obtained bv lDuez et al. I (1201 Oh . in which Newtonian 
resistive MHD simulations are performed. These results seem 
to be incompatible with our present results . The reason f or this 
disc repancy seems to be the following: Lan der et al. I (1201 Oh 
and ILander & Jones ] (12010) restrict their studies to low-order 
azimuthal modes (m < 6) and lDuez et al. I (1201 Oh employs some 
kinds of (artificial) viscosity and resistivity for the evolution to 
remove numerical instability caused by a shor t-wavelength os- 
cillation. Note that Lander et al. I (12010) and Lander & Jones I 
(i201Q) also use artificial viscosity to stabilize their computa- 
tions. This suggests that in their simulations, short-wavelength 
modes might be suppressed. As can be seen from Figures [T] 
and [3] the Parker instability found in this study is characterized 
by a high-order azimuthal mode number, which implies that the 
unstable modes of t he Parker instabihty near the surface have a 
short wavelength. In lDuez et al. I (1201 Oh . they used a stably strat- 
ified model with the polytrope index n - 3 and forecasted the 
instability as we did in Figure [T] They found the primarily in- 
stability is not the Parker instability, but the Taylor instability in 
their model. This is likely that the restoring force due to the en- 
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tropy gradient near the surface stabilizes the magnetic buoyant 
force. 

Finally, we comment on a possible potential effect for stabi- 
lizing the magnetic field instabilities which are not taken into ac- 
count in our present work. In a stably stratified region of the star, 
the buoyant force inhibits the interchange , Parker, and Tayler 
instabilities from growing (Ach eson"lll978 ). For the cold neu- 
tron star, the composition gradient induced by chemical inhomo- 
geneities can stably stratify the neutron star matter a nd supports 
the gravity mode (g-mode) oscillations. iFinn 1 (fT987l) considered 
crustal g-modes due to the composition discontinuities in the 
outer envelopes, whose typical oscillation per iod is about 5 ms 
for a canonical neutron star model. Reis enegger & Goldreich I 
(Il992h studied the effect of g-modes associated with buoyancy 
induced by proton-neuron composition gradients in the core, 
whose typical oscillation period is about 2 ms for a canonical 
neutron star model. If the growth time scale of the Parker in- 
stability found in our study is longer than these periods of the 
g-modes, it is possible that the buoyant forces inside the neu- 
tron star suppress the growth of the Parker instability. In the 
present study, we find that a typical growth time scale of the 
Parker instability is about 20Mo, which gives a typical growth 
time scale of 0.1 ms for a canonical neutron star model with a 
strong magnetic field of 10'^ G. For the present models, thus, it 
seems that the buoyancy inside the neutron star could not sup- 
press the onset of the Parker instability. However, for a weaker 
magnetic field strength » 10''' G, the Parker instability could be 
suppressed by the buoyancy because a typical growth time scale 
of the Parker instability becomes 1 ms. To derive a definite an- 
swer to this problem, whether or not the buoyancy can stabilize 
the Parker instability inside the neutron star, we have to know 
the local growth rate of the Parker instability. Unfortunately, in 
our simulations, it is difficult to estimate it in the vicinity of the 
stellar surface because the local magnetic structure highly de- 
pends on the numerical resolution. To investigate this point more 
precisely, we have to take into account the chemical inhomo- 
geneities; this implies that it is necessary to implement equations 
of state which depend on the chemical compositions and to ob- 
tain the evolution of the chemical compositions. This is beyond 
the scope of this paper. 
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Fig. 1. The growth rate of the instability of a neutron star with toroidal magnetic fields normalized by the maximum growth rate 
Icmaxl (left) and corresponding m mode (right) on the meridian plane for non-rotating Newtonian polytrope with po = 0.16 and 
Bo - 0.2. The white colored region and the dotted curve denote the stable region and critical radius, respectively, rj denotes the 
stellar radius. 
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Fig. 11. The same as Figure[TO]but for model R22H2T8. 
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